How does microbial community structure change with depth and oxygen concentration?

Alpha-diversity

Alpha-diversity and richness were calculated for the total community in R.

load("qiime2_phyloseq.RData")

Samples were rarefied/normalized to 100,000 sequences per sample to facilitate comparisons between samples. A random seed was set to ensure reproducibility.

set.seed(4832)
m.norm = rarefy_even_depth(qiime2, sample.size=100000)
## You set `rngseed` to FALSE. Make sure you've set & recorded
##  the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## 6OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...

Rarefied counts were converted to relative abundance percentages.

m.perc = transform_sample_counts(m.norm, function(x) 100 * x/sum(x))
# Calculation
m.alpha = estimate_richness(m.norm, measures = c("Chao1", "Shannon"))
m.meta.alpha = full_join(rownames_to_column(m.alpha), rownames_to_column(data.frame(m.perc@sam_data)), by = "rowname")
m.meta.alpha

Example plots

m.meta.alpha %>% 
  
ggplot() +
  geom_point(aes(x=Depth_m, y=Shannon)) +
   geom_smooth(method='auto', aes(x=as.numeric(Depth_m), y=Shannon)) +
  labs(title="Alpha-diversity across depth", y="Shannon's diversity index", x="Depth (m)")
## `geom_smooth()` using method = 'loess'

m.meta.alpha %>% 

  ggplot() +
  geom_point(aes(x=O2_uM, y=Shannon)) +
   geom_smooth(method='auto', aes(x=as.numeric(O2_uM), y=Shannon)) +
  labs(title="Alpha-diversity across oxygen", y="Shannon's diversity index", x="O2 (uM)")
## `geom_smooth()` using method = 'loess'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -1.0833
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 33.437
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 0
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1046.8
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at -1.0833
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 33.437
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 0
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 1046.8

m.meta.alpha %>%
  
  ggplot() +
  geom_point(aes(x=O2_uM, y=Shannon)) +
  labs(title="Alpha-diversity across oxygen", y="Shannon's diversity index", x="Oxygen (uM)")

m.meta.alpha %>% 
  mutate(O2_group = ifelse(O2_uM == 0, "anoxic", "oxic")) %>% 

# Use this in a plot
ggplot() +
  geom_boxplot(aes(x=O2_group, y=Shannon)) +
  labs(title="Alpha-diversity by oxic/anoxic", y="Shannon's diversity index", x="Oxygen")

Taxa presence and abundance

Example plots using phyloseq for domain data. You should explore other taxonomy levels.

m.perc %>% 
plot_bar(fill="Domain") + 
  geom_bar(aes(fill=Domain), stat="identity") +
  labs(title="Domains across samples") +
  labs(fill='Domain') +
  theme(plot.title = element_text(size=35, hjust=.5,vjust=.5,face="plain"),
        strip.text = element_text(size=12),
        axis.title.x = element_text(size=25),
        axis.title.y = element_text(size=25),
        axis.text.x = element_text(colour="grey20",size=12,angle=45,hjust=1,vjust=1),
        axis.text.y = element_text(colour="grey20",size=12),
        legend.text = element_text(size=12),
        legend.box.background = element_rect(),
        legend.box.margin = margin(6, 6, 6, 6),
        legend.position = "right")

m.perc %>% 
plot_bar() + 
  geom_bar(aes(fill=Domain), stat="identity") +
  facet_wrap(~Phylum, scales="free_y")+
  labs(title="Phyla across samples")+
  labs(fill='Domain') +
  theme(plot.title = element_text(size=35, hjust=.5,vjust=.5,face="plain"),
        strip.text = element_text(size=5),
        axis.title.x = element_text(size=25),
        axis.title.y = element_text(size=25),
        axis.text.x = element_text(colour="grey20",size=9,angle=45,hjust=1,vjust=1),
        axis.text.y = element_text(colour="grey20",size=9),
        legend.text = element_text(size=12),
        legend.box.background = element_rect(),
        legend.box.margin = margin(6, 6, 6, 6),
        legend.position = "bottom")

Example plots using phyloseq for Family data.

# load subset
m.perc %>% 
plot_bar(fill="Family") + 
  geom_bar(aes(fill=Family), stat="identity") +
  labs(title="Family across samples") +
  theme(plot.title = element_text(size=17.5, hjust=.5,vjust=.5,face="plain"),
        axis.title.x = element_text(size=12.5),
        axis.title.y = element_text(size=12.5),
        axis.text.x = element_text(colour="grey20",size=7.5,angle=45,hjust=1,vjust=1),
        axis.text.y = element_text(colour="grey20",size=7.5),
        legend.text = element_text(size=6),
        legend.box.background = element_rect(),
        legend.box.margin = margin(6, 6, 6, 6),
        legend.position = "bottom")

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
plot_bar(fill="Genus") + 
  geom_bar(aes(fill=Genus), stat="identity") +
  labs(title="Genus of family Oceanospirillaceae across samples") +
  labs(fill="Genus") +
  theme(plot.title = element_text(size=17.5, hjust=.5,vjust=.5,face="plain"),
        axis.title.x = element_text(size=12.5),
        axis.title.y = element_text(size=12.5),
        axis.text.x = element_text(colour="grey20",size=7.5,angle=45,hjust=1,vjust=1),
        axis.text.y = element_text(colour="grey20",size=7.5),
        legend.text = element_text(size=6),
        legend.box.background = element_rect(),
        legend.box.margin = margin(6, 6, 6, 6),
        legend.position = "bottom")

Examples outside of phyloseq

m.perc %>%
  tax_glom(taxrank = 'Domain') %>%
  psmelt() %>% 

ggplot() +
  geom_boxplot(aes(x=Domain, y=Abundance)) +
  coord_flip() +
  labs(title="Domain boxplots")+
  theme(plot.title = element_text(size=17.5, hjust=.5,vjust=.5,face="plain"),
      axis.title.x = element_text(size=12.5),
      axis.title.y = element_text(size=12.5)
  )

m.perc %>% 
  tax_glom(taxrank = 'Family') %>%
  psmelt() %>% 

ggplot() +
  geom_boxplot(aes(x=Family, y=Abundance)) +
  coord_flip() +
  labs(title="Family boxplots")+
  theme(plot.title = element_text(size=17.5, hjust=.5,vjust=.5,face="plain"),
    axis.title.x = element_text(size=12.5),
    axis.title.y = element_text(size=12.5)
  )

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
  psmelt() %>% 

ggplot() +
  geom_boxplot(aes(x=Genus, y=Abundance)) +
  coord_flip() +
  labs(title="Domain boxplots")+
  theme(plot.title = element_text(size=17.5, hjust=.5,vjust=.5,face="plain"),
    axis.title.x = element_text(size=12.5),
    axis.title.y = element_text(size=12.5)
  )

Does your taxon of interest significantly differ in abundance with depth and/or oxygen concentration?

Using the magrittr package, we can pipe our tidyverse modified data into linear models and other statiscal tests.

Depth Linear model

m.norm %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
  tax_glom(taxrank = 'Family') %>%
  psmelt() %>%

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       6       5       7 
## -152.74  545.73   40.05 -200.21   42.28 -344.46   69.35 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2405.898    290.936   8.270 0.000422 ***
## Depth_m      -11.516      2.115  -5.445 0.002838 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 312.4 on 5 degrees of freedom
## Multiple R-squared:  0.8557, Adjusted R-squared:  0.8268 
## F-statistic: 29.65 on 1 and 5 DF,  p-value: 0.002838

Plot to go along with linear model above.

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
  psmelt() %>% 
  group_by(Sample) %>% 
  summarize(Abundance_sum=sum(Abundance), Depth_m=mean(Depth_m)) %>% 
  
ggplot() +
  geom_point(aes(x=Depth_m, y=Abundance_sum)) +
  geom_smooth(method='lm', aes(x=as.numeric(Depth_m), y=Abundance_sum)) +
  labs(title="Abundance of Family Oceanospirillaceae across depth")

Oxygen Linear model

m.norm %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
  tax_glom(taxrank = 'Family') %>%
  psmelt() %>%

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       6       5       7 
## -173.54  888.19  196.52 -123.18  -65.99 -279.99 -441.99 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  613.993    212.371   2.891   0.0341 *
## O2_uM          7.835      2.516   3.113   0.0264 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 479.8 on 5 degrees of freedom
## Multiple R-squared:  0.6597, Adjusted R-squared:  0.5917 
## F-statistic: 9.693 on 1 and 5 DF,  p-value: 0.02645

Plot to go along with linear model above.

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
  psmelt() %>% 
  group_by(Sample) %>% 
  summarize(Abundance_sum=sum(Abundance), O2_uM=mean(O2_uM)) %>% 
  
ggplot() +
  geom_point(aes(x=O2_uM, y=Abundance_sum)) +
  geom_smooth(method='lm', aes(x=as.numeric(O2_uM), y=Abundance_sum)) +
  labs(title="Abundance of Family Oceanospirillaceaes across oxygen")

Within your taxon, what is the richness (number of OTUs/ASVs)?

Across all samples, there are 28 OTUs (i.e. taxa) within family Oceanospirillaceae

m.norm %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 28 taxa and 7 samples ]
## sample_data() Sample Data:       [ 7 samples by 22 sample variables ]
## tax_table()   Taxonomy Table:    [ 28 taxa by 7 taxonomic ranks ]

Counting OTUs within family Oceanospirillaceae within each sample…

m.norm %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>%
  estimate_richness(measures = c("Observed"))
## Warning in estimate_richness(., measures = c("Observed")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.

Do the abundances of OTUs/ASVs within your taxon of interest change significantly with depth and/or oxygen concentration?

Depth General linear model for each OTU

#Asv107
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv107") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  4.7185 -0.7954 -1.1195 -0.9574  0.3928 -0.5794 -1.6596 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.50082    2.21579  -0.226    0.830
## Depth_m      0.01080    0.01611   0.671    0.532
## 
## Residual standard error: 2.38 on 5 degrees of freedom
## Multiple R-squared:  0.08252,    Adjusted R-squared:  -0.101 
## F-statistic: 0.4497 on 1 and 5 DF,  p-value: 0.5322
#Asv120
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv120") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  25.145  25.560   1.124 -14.617 -10.481 -13.202 -13.529 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.26350   18.47249   0.556    0.602
## Depth_m      0.02177    0.13429   0.162    0.878
## 
## Residual standard error: 19.84 on 5 degrees of freedom
## Multiple R-squared:  0.005227,   Adjusted R-squared:  -0.1937 
## F-statistic: 0.02627 on 1 and 5 DF,  p-value: 0.8776
#Asv131
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv131") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  10.2285   2.7496  10.5447  -0.5912 -11.7270  -7.2727  -3.9319 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 33.99869    8.63139   3.939   0.0110 *
## Depth_m     -0.22272    0.06275  -3.549   0.0164 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.269 on 5 degrees of freedom
## Multiple R-squared:  0.7159, Adjusted R-squared:  0.6591 
## F-statistic:  12.6 on 1 and 5 DF,  p-value: 0.0164
#Asv213
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv213") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  215.21 -238.29 -161.73  -82.81   73.02  -25.89  220.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  845.073    181.870   4.647   0.0056 **
## Depth_m       -5.328      1.322  -4.030   0.0100 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 195.3 on 5 degrees of freedom
## Multiple R-squared:  0.7646, Adjusted R-squared:  0.7175 
## F-statistic: 16.24 on 1 and 5 DF,  p-value: 0.01002
#Asv243
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv243") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## 135.300  48.178   5.336 -16.347 -40.505 -52.311 -79.651 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  81.0897    74.0794   1.095    0.324
## Depth_m      -0.1439     0.5385  -0.267    0.800
## 
## Residual standard error: 79.55 on 5 degrees of freedom
## Multiple R-squared:  0.01408,    Adjusted R-squared:  -0.1831 
## F-statistic: 0.07139 on 1 and 5 DF,  p-value: 0.8
#Asv417
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv417") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 17.113 -2.651 -3.463 -2.468 -2.809 -2.992 -2.730 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  3.515548   7.704200   0.456    0.667
## Depth_m     -0.005237   0.056008  -0.094    0.929
## 
## Residual standard error: 8.274 on 5 degrees of freedom
## Multiple R-squared:  0.001746,   Adjusted R-squared:  -0.1979 
## F-statistic: 0.008744 on 1 and 5 DF,  p-value: 0.9291
#Asv476
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv476") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  24.170  -4.936  -6.979 -11.064  -2.213  -9.021  10.043 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.40426   12.94853  -0.881    0.419
## Depth_m       0.13617    0.09413   1.447    0.208
## 
## Residual standard error: 13.91 on 5 degrees of freedom
## Multiple R-squared:  0.295,  Adjusted R-squared:  0.154 
## F-statistic: 2.093 on 1 and 5 DF,  p-value: 0.2077
#Asv668
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv668") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 15.106 -1.383 -6.915 -5.638 -4.362 -3.085  6.277 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.12766    8.09283  -0.881    0.419
## Depth_m      0.08511    0.05883   1.447    0.208
## 
## Residual standard error: 8.691 on 5 degrees of freedom
## Multiple R-squared:  0.295,  Adjusted R-squared:  0.154 
## F-statistic: 2.093 on 1 and 5 DF,  p-value: 0.2077
#Asv810
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv810") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  158.10 -165.76  -42.37 -131.39  -23.34  177.07   27.68 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 676.5840   134.4226   5.033  0.00399 **
## Depth_m      -4.2683     0.9772  -4.368  0.00724 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 144.4 on 5 degrees of freedom
## Multiple R-squared:  0.7923, Adjusted R-squared:  0.7508 
## F-statistic: 19.08 on 1 and 5 DF,  p-value: 0.007237
#Asv1050
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1050") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 13.690 -2.247 -2.121 -2.771 -2.184 -2.393 -1.974 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.81244    6.16336   0.456    0.667
## Depth_m     -0.00419    0.04481  -0.094    0.929
## 
## Residual standard error: 6.619 on 5 degrees of freedom
## Multiple R-squared:  0.001746,   Adjusted R-squared:  -0.1979 
## F-statistic: 0.008744 on 1 and 5 DF,  p-value: 0.9291
#Asv1152
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1152") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  68.962   5.881  -1.242 -29.669 -11.447 -14.804 -17.681 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.62848   33.14541   0.924    0.398
## Depth_m     -0.09591    0.24096  -0.398    0.707
## 
## Residual standard error: 35.59 on 5 degrees of freedom
## Multiple R-squared:  0.03071,    Adjusted R-squared:  -0.1631 
## F-statistic: 0.1584 on 1 and 5 DF,  p-value: 0.707
#Asv1283
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1283") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  69.31  53.87  37.03 -15.12 -17.55 -54.57 -72.97 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 107.2812    55.7393   1.925    0.112
## Depth_m      -0.1716     0.4052  -0.423    0.690
## 
## Residual standard error: 59.86 on 5 degrees of freedom
## Multiple R-squared:  0.03461,    Adjusted R-squared:  -0.1585 
## F-statistic: 0.1792 on 1 and 5 DF,  p-value: 0.6896
#Asv1352
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1352") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 152.08  95.68 -24.61 -23.21 -39.91 -62.90 -97.13 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  98.9349    91.6288    1.08    0.330
## Depth_m      -0.1802     0.6661   -0.27    0.798
## 
## Residual standard error: 98.4 on 5 degrees of freedom
## Multiple R-squared:  0.01442,    Adjusted R-squared:  -0.1827 
## F-statistic: 0.07315 on 1 and 5 DF,  p-value: 0.7976
#Asv1414
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1414") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -42.235  83.392  55.636   1.532 -45.259 -51.155  -1.911 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 283.3051    53.5000   5.295   0.0032 **
## Depth_m      -1.4070     0.3889  -3.618   0.0153 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 57.45 on 5 degrees of freedom
## Multiple R-squared:  0.7235, Adjusted R-squared:  0.6683 
## F-statistic: 13.09 on 1 and 5 DF,  p-value: 0.01526
#Asv1430
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1430") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  23.566   9.791  -4.813  -2.157  -6.804  -8.796 -10.787 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.11915   12.62482  -0.881    0.419
## Depth_m       0.13277    0.09178   1.447    0.208
## 
## Residual standard error: 13.56 on 5 degrees of freedom
## Multiple R-squared:  0.295,  Adjusted R-squared:  0.154 
## F-statistic: 2.093 on 1 and 5 DF,  p-value: 0.2077
#Asv1486
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1486") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  38.201  10.489  -5.357  -7.612  -8.578 -17.597  -9.545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.24124   19.22078   0.949    0.386
## Depth_m     -0.06442    0.13973  -0.461    0.664
## 
## Residual standard error: 20.64 on 5 degrees of freedom
## Multiple R-squared:  0.04077,    Adjusted R-squared:  -0.1511 
## F-statistic: 0.2125 on 1 and 5 DF,  p-value: 0.6641
#Asv1490
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1490") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  43.642   9.631 -10.287  -8.298  -9.828 -10.746 -14.113 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.41899   21.16585   0.681    0.526
## Depth_m     -0.03061    0.15387  -0.199    0.850
## 
## Residual standard error: 22.73 on 5 degrees of freedom
## Multiple R-squared:  0.00785,    Adjusted R-squared:  -0.1906 
## F-statistic: 0.03956 on 1 and 5 DF,  p-value: 0.8502
#Asv1570
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1570") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  78.947  13.491  10.544   2.439 -13.613 -48.740 -43.068 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  69.0386    43.5720   1.584    0.174
## Depth_m      -0.1299     0.3168  -0.410    0.699
## 
## Residual standard error: 46.79 on 5 degrees of freedom
## Multiple R-squared:  0.03252,    Adjusted R-squared:  -0.161 
## F-statistic: 0.1681 on 1 and 5 DF,  p-value: 0.6988
#Asv1593
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1593") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 32.515 -5.187 -5.684 -6.580 -5.038 -5.336 -4.689 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  6.679542  14.637980   0.456    0.667
## Depth_m     -0.009951   0.106415  -0.094    0.929
## 
## Residual standard error: 15.72 on 5 degrees of freedom
## Multiple R-squared:  0.001746,   Adjusted R-squared:  -0.1979 
## F-statistic: 0.008744 on 1 and 5 DF,  p-value: 0.9291
#Asv1797
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1797") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 19.336 -8.851 -5.583 -7.217 -1.770 -3.949  8.034 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.12340   10.35883  -0.881    0.419
## Depth_m      0.10894    0.07531   1.447    0.208
## 
## Residual standard error: 11.12 on 5 degrees of freedom
## Multiple R-squared:  0.295,  Adjusted R-squared:  0.154 
## F-statistic: 2.093 on 1 and 5 DF,  p-value: 0.2077
#Asv1977
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1977") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 14.502 -4.187  6.026 -1.328 -6.638 -5.413 -2.962 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.84255    7.76912  -0.881    0.419
## Depth_m      0.08170    0.05648   1.447    0.208
## 
## Residual standard error: 8.343 on 5 degrees of freedom
## Multiple R-squared:  0.295,  Adjusted R-squared:  0.154 
## F-statistic: 2.093 on 1 and 5 DF,  p-value: 0.2077
#Asv2036
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2036") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 26.449 -4.947 -2.902 -4.749 -4.353 -4.089 -5.409 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.77021   11.92437   0.232    0.826
## Depth_m      0.01319    0.08669   0.152    0.885
## 
## Residual standard error: 12.81 on 5 degrees of freedom
## Multiple R-squared:  0.00461,    Adjusted R-squared:  -0.1945 
## F-statistic: 0.02316 on 1 and 5 DF,  p-value: 0.885
#Asv2097
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2097") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 109.23  62.72 -12.41 -21.55 -21.68 -69.98 -46.32 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  71.2291    64.5079   1.104    0.320
## Depth_m      -0.1245     0.4690  -0.266    0.801
## 
## Residual standard error: 69.28 on 5 degrees of freedom
## Multiple R-squared:  0.01391,    Adjusted R-squared:  -0.1833 
## F-statistic: 0.07054 on 1 and 5 DF,  p-value: 0.8012
#Asv2125
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2125") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  52.918 -17.870  -6.180  -7.311  -8.443  -9.574  -3.540 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.62390   24.23491   0.768    0.477
## Depth_m     -0.07542    0.17618  -0.428    0.686
## 
## Residual standard error: 26.03 on 5 degrees of freedom
## Multiple R-squared:  0.03535,    Adjusted R-squared:  -0.1576 
## F-statistic: 0.1832 on 1 and 5 DF,  p-value: 0.6864
#Asv2215
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2215") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  9.0953 -1.4511 -1.0622 -0.6085 -3.0714 -1.2566 -1.6455 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  3.20098    4.16538   0.768    0.477
## Depth_m     -0.01296    0.03028  -0.428    0.686
## 
## Residual standard error: 4.473 on 5 degrees of freedom
## Multiple R-squared:  0.03535,    Adjusted R-squared:  -0.1576 
## F-statistic: 0.1832 on 1 and 5 DF,  p-value: 0.6864
#Asv2281
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2281") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## 102.405  70.715   9.423  -4.813 -23.049 -63.179 -91.502 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  61.6887    70.2964   0.878    0.420
## Depth_m       0.1491     0.5110   0.292    0.782
## 
## Residual standard error: 75.49 on 5 degrees of freedom
## Multiple R-squared:  0.01673,    Adjusted R-squared:  -0.1799 
## F-statistic: 0.08508 on 1 and 5 DF,  p-value: 0.7822
#Asv2297
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2297") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  33.147  10.588   9.581 -15.839 -13.081 -12.750 -11.647 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.05925   18.72956   0.857    0.430
## Depth_m     -0.02206    0.13616  -0.162    0.878
## 
## Residual standard error: 20.11 on 5 degrees of freedom
## Multiple R-squared:  0.005223,   Adjusted R-squared:  -0.1937 
## F-statistic: 0.02625 on 1 and 5 DF,  p-value: 0.8776
#Asv2341
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2341") %>% 

  lm(Abundance ~ Depth_m, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ Depth_m, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  7.2511 -2.0936 -1.4809 -3.3191 -2.7064  3.0128 -0.6638 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.42128    3.88456  -0.881    0.419
## Depth_m      0.04085    0.02824   1.447    0.208
## 
## Residual standard error: 4.172 on 5 degrees of freedom
## Multiple R-squared:  0.295,  Adjusted R-squared:  0.154 
## F-statistic: 2.093 on 1 and 5 DF,  p-value: 0.2077

Repeated for all OTUs within domain and then corrected for multiple comparisons.

p.adjust(c(0.5322, 0.8776, 0.0164, 0.01002, 0.8, 0.9291, 0.2077, 0.2077, 0.007237, 0.9291, 0.707, 0.6896, 0.7976, 0.01526, 0.2077, 0.6641, 0.8502, 0.6988, 0.9291, 0.2077, 0.2077, 0.885, 0.8012, 0.6864, 0.6864, 0.7822, 0.8776, 0.2077), method="fdr")
##  [1] 0.92910 0.92910 0.11480 0.11480 0.92910 0.92910 0.58156 0.58156
##  [9] 0.11480 0.92910 0.92910 0.92910 0.92910 0.11480 0.58156 0.92910
## [17] 0.92910 0.92910 0.92910 0.58156 0.58156 0.92910 0.92910 0.92910
## [25] 0.92910 0.92910 0.92910 0.58156

Plots to go along with these tests.

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
  psmelt() %>% 
  
ggplot() +
  geom_point(aes(x=Depth_m, y=Abundance)) +
  geom_smooth(method='lm', aes(x=Depth_m, y=Abundance)) +
  facet_wrap(~OTU, scales="free_y") +
  labs(title="Example 11: Abundance of OTUs within Family Oceanospirillaceae across depth")

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>%
  psmelt() %>% 
  
ggplot() +
  geom_point(aes(x=Sample, y=OTU, size=Abundance, color=OTU)) + 
  scale_size_continuous(range = c(0,5)) +
  labs(y="ASV", title="Abundance of ASVs within Family Oceanospirillaceae across depth")+
  theme(plot.title = element_text(size=15, hjust=.5,vjust=.5,face="plain"),
        axis.text.x = element_text(colour="grey20",size=10,angle=45,hjust=.5,vjust=.5,face="plain"),
        legend.box.background = element_rect(),
        legend.box.margin = margin(6, 6, 6, 6),
        legend.position = "right")

Oxygen

#Asv107
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv107") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  4.8244 -0.9410 -1.1756 -1.0274  0.3954 -0.9000 -1.1756 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.175639   1.065055   1.104     0.32
## O2_uM       -0.007251   0.012620  -0.575     0.59
## 
## Residual standard error: 2.406 on 5 degrees of freedom
## Multiple R-squared:  0.06193,    Adjusted R-squared:  -0.1257 
## F-statistic: 0.3301 on 1 and 5 DF,  p-value: 0.5905
#Asv120
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv120") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##        1        2        3        4        5        6        7 
##  23.4624  24.6584   0.3315 -15.5376  -3.0208 -14.3564 -15.5376 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.53756    8.52995   1.822    0.128
## O2_uM       -0.05777    0.10108  -0.572    0.592
## 
## Residual standard error: 19.27 on 5 degrees of freedom
## Multiple R-squared:  0.06133,    Adjusted R-squared:  -0.1264 
## F-statistic: 0.3267 on 1 and 5 DF,  p-value: 0.5924
#Asv131
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv131") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  1.521  2.768  2.768  2.768 -4.820 -3.690 -1.313 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.76761    1.58199  -1.749 0.140622    
## O2_uM        0.19960    0.01875  10.648 0.000126 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.574 on 5 degrees of freedom
## Multiple R-squared:  0.9578, Adjusted R-squared:  0.9493 
## F-statistic: 113.4 on 1 and 5 DF,  p-value: 0.0001264
#Asv213
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv213") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  26.64 -73.73 -77.36 -22.86  68.43  49.43  29.43 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -29.432     28.443  -1.035    0.348    
## O2_uM          4.661      0.337  13.828 3.55e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 64.26 on 5 degrees of freedom
## Multiple R-squared:  0.9745, Adjusted R-squared:  0.9694 
## F-statistic: 191.2 on 1 and 5 DF,  p-value: 3.552e-05
#Asv243
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv243") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##        1        2        3        4        5        6        7 
## 137.9101  46.8673  -0.3275 -30.0958 -52.0958 -71.0958 -31.1625 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  71.0958    34.7763   2.044   0.0963 .
## O2_uM        -0.1843     0.4121  -0.447   0.6734  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 78.56 on 5 degrees of freedom
## Multiple R-squared:  0.03847,    Adjusted R-squared:  -0.1538 
## F-statistic:   0.2 on 1 and 5 DF,  p-value: 0.6734
#Asv417
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv417") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 17.069 -3.137 -1.757 -3.137 -3.007 -2.895 -3.137 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  3.136820   3.657612   0.858    0.430
## O2_uM       -0.006367   0.043341  -0.147    0.889
## 
## Residual standard error: 8.263 on 5 degrees of freedom
## Multiple R-squared:  0.004298,   Adjusted R-squared:  -0.1948 
## F-statistic: 0.02158 on 1 and 5 DF,  p-value: 0.8889
#Asv476
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv476") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 32.162 -6.274 -6.849 -7.838 -6.000 -7.838  2.636 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  7.83759    7.10036   1.104     0.32
## O2_uM       -0.04834    0.08414  -0.575     0.59
## 
## Residual standard error: 16.04 on 5 degrees of freedom
## Multiple R-squared:  0.06193,    Adjusted R-squared:  -0.1257 
## F-statistic: 0.3301 on 1 and 5 DF,  p-value: 0.5905
#Asv668
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv668") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 20.102 -3.750 -4.898 -4.898 -4.281 -3.921  1.647 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  4.89850    4.43773   1.104     0.32
## O2_uM       -0.03021    0.05258  -0.575     0.59
## 
## Residual standard error: 10.03 on 5 degrees of freedom
## Multiple R-squared:  0.06193,    Adjusted R-squared:  -0.1257 
## F-statistic: 0.3301 on 1 and 5 DF,  p-value: 0.5905
#Asv810
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv810") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  15.177 -34.199   4.559 -64.341  34.934  21.934  21.934 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -21.9345    17.4194  -1.259    0.264    
## O2_uM         3.6866     0.2064  17.860 1.01e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.35 on 5 degrees of freedom
## Multiple R-squared:  0.9846, Adjusted R-squared:  0.9815 
## F-statistic:   319 on 1 and 5 DF,  p-value: 1.01e-05
#Asv1050
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1050") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 13.655 -2.405 -2.509 -1.406 -2.509 -2.316 -2.509 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.509456   2.926090   0.858    0.430
## O2_uM       -0.005094   0.034673  -0.147    0.889
## 
## Residual standard error: 6.61 on 5 degrees of freedom
## Multiple R-squared:  0.004298,   Adjusted R-squared:  -0.1948 
## F-statistic: 0.02158 on 1 and 5 DF,  p-value: 0.8889
#Asv1152
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1152") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  71.188   5.957  -5.360 -11.536 -20.360 -20.360 -19.528 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.36044   15.92980   1.278    0.257
## O2_uM       -0.04073    0.18876  -0.216    0.838
## 
## Residual standard error: 35.99 on 5 degrees of freedom
## Multiple R-squared:  0.009225,   Adjusted R-squared:  -0.1889 
## F-statistic: 0.04656 on 1 and 5 DF,  p-value: 0.8377
#Asv1283
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1283") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  69.21  57.74  26.21 -18.89 -25.79 -18.68 -89.79 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 89.79117   26.73869   3.358   0.0201 *
## O2_uM       -0.09281    0.31684  -0.293   0.7813  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 60.41 on 5 degrees of freedom
## Multiple R-squared:  0.01687,    Adjusted R-squared:  -0.1798 
## F-statistic: 0.08581 on 1 and 5 DF,  p-value: 0.7813
#Asv1352
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1352") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 155.42  94.17 -31.44 -39.93 -53.93 -85.93 -38.36 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  85.9292    43.0850   1.994    0.103
## O2_uM        -0.2195     0.5105  -0.430    0.685
## 
## Residual standard error: 97.33 on 5 degrees of freedom
## Multiple R-squared:  0.03567,    Adjusted R-squared:  -0.1572 
## F-statistic: 0.1849 on 1 and 5 DF,  p-value: 0.6851
#Asv1414
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1414") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## -30.79 124.75  63.14  19.71 -40.94 -67.94 -67.94 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  67.9400    35.2983   1.925   0.1122  
## O2_uM         0.8762     0.4183   2.095   0.0903 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 79.74 on 5 degrees of freedom
## Multiple R-squared:  0.4674, Adjusted R-squared:  0.3609 
## F-statistic: 4.389 on 1 and 5 DF,  p-value: 0.09034
#Asv1430
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1430") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 31.358  2.570 -6.117 -5.850 -6.678 -7.642 -7.642 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  7.64165    6.92285   1.104     0.32
## O2_uM       -0.04713    0.08203  -0.575     0.59
## 
## Residual standard error: 15.64 on 5 degrees of freedom
## Multiple R-squared:  0.06193,    Adjusted R-squared:  -0.1257 
## F-statistic: 0.3301 on 1 and 5 DF,  p-value: 0.5905
#Asv1486
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1486") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  39.770  10.686 -10.794 -10.794 -10.794  -7.583 -10.491 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.79384    9.31225   1.159    0.299
## O2_uM       -0.01482    0.11034  -0.134    0.898
## 
## Residual standard error: 21.04 on 5 degrees of freedom
## Multiple R-squared:  0.003595,   Adjusted R-squared:  -0.1957 
## F-statistic: 0.01804 on 1 and 5 DF,  p-value: 0.8984
#Asv1490
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1490") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  44.240   7.027 -11.321 -11.973 -11.973 -10.941  -5.059 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.97305   10.02915   1.194    0.286
## O2_uM       -0.03191    0.11884  -0.269    0.799
## 
## Residual standard error: 22.66 on 5 degrees of freedom
## Multiple R-squared:  0.01421,    Adjusted R-squared:  -0.1829 
## F-statistic: 0.07209 on 1 and 5 DF,  p-value: 0.799
#Asv1570
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1570") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  81.706   9.986  10.152  -5.017 -23.017 -16.793 -57.017 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 57.01715   20.73066   2.750   0.0403 *
## O2_uM       -0.09796    0.24565  -0.399   0.7065  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 46.83 on 5 degrees of freedom
## Multiple R-squared:  0.03082,    Adjusted R-squared:  -0.163 
## F-statistic: 0.159 on 1 and 5 DF,  p-value: 0.7065
#Asv1593
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1593") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 32.431 -5.960 -5.500 -3.339 -5.960 -5.713 -5.960 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  5.95996    6.94946   0.858    0.430
## O2_uM       -0.01210    0.08235  -0.147    0.889
## 
## Residual standard error: 15.7 on 5 degrees of freedom
## Multiple R-squared:  0.004298,   Adjusted R-squared:  -0.1948 
## F-statistic: 0.02158 on 1 and 5 DF,  p-value: 0.8889
#Asv1797
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1797") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 25.730 -6.270 -5.479 -6.270 -4.800 -5.019  2.109 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  6.27007    5.68029   1.104     0.32
## O2_uM       -0.03867    0.06731  -0.575     0.59
## 
## Residual standard error: 12.83 on 5 degrees of freedom
## Multiple R-squared:  0.06193,    Adjusted R-squared:  -0.1257 
## F-statistic: 0.3301 on 1 and 5 DF,  p-value: 0.5905
#Asv1977
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv1977") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 19.297 -4.110  1.581 -3.600 -4.703 -4.703 -3.764 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  4.70256    4.26022   1.104     0.32
## O2_uM       -0.02900    0.05048  -0.575     0.59
## 
## Residual standard error: 9.624 on 5 degrees of freedom
## Multiple R-squared:  0.06193,    Adjusted R-squared:  -0.1257 
## F-statistic: 0.3301 on 1 and 5 DF,  p-value: 0.5905
#Asv2036
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2036") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
## 26.1013 -5.3082 -0.9694 -5.3082 -4.6603 -4.5470 -5.3082 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  5.30818    5.63103   0.943    0.389
## O2_uM       -0.02002    0.06672  -0.300    0.776
## 
## Residual standard error: 12.72 on 5 degrees of freedom
## Multiple R-squared:  0.01769,    Adjusted R-squared:  -0.1788 
## F-statistic: 0.09007 on 1 and 5 DF,  p-value: 0.7762
#Asv2097
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2097") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 111.45  61.52 -17.44 -31.81 -33.81 -27.09 -62.81 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  62.8128    30.2469   2.077   0.0925 .
## O2_uM        -0.1649     0.3584  -0.460   0.6649  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 68.33 on 5 degrees of freedom
## Multiple R-squared:  0.0406, Adjusted R-squared:  -0.1513 
## F-statistic: 0.2116 on 1 and 5 DF,  p-value: 0.6649
#Asv2125
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2125") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
## 54.796 -7.344 -9.600 -9.600 -9.387 -9.263 -9.600 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  9.60022   11.72301   0.819    0.450
## O2_uM       -0.01041    0.13891  -0.075    0.943
## 
## Residual standard error: 26.48 on 5 degrees of freedom
## Multiple R-squared:  0.001122,   Adjusted R-squared:  -0.1987 
## F-statistic: 0.005619 on 1 and 5 DF,  p-value: 0.9432
#Asv2215
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2215") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  9.418 -1.613 -1.650 -1.650 -1.262 -1.650 -1.592 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.65004    2.01489   0.819    0.450
## O2_uM       -0.00179    0.02388  -0.075    0.943
## 
## Residual standard error: 4.552 on 5 degrees of freedom
## Multiple R-squared:  0.001122,   Adjusted R-squared:  -0.1987 
## F-statistic: 0.005619 on 1 and 5 DF,  p-value: 0.9432
#Asv2281
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2281") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  96.380  60.291   4.283 -12.131 -35.709 -16.404 -96.709 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  96.7092    30.6591   3.154   0.0253 *
## O2_uM        -0.3706     0.3633  -1.020   0.3544  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 69.26 on 5 degrees of freedom
## Multiple R-squared:  0.1723, Adjusted R-squared:  0.006764 
## F-statistic: 1.041 on 1 and 5 DF,  p-value: 0.3544
#Asv2297
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2297") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  33.467  10.230   6.875  -6.052 -14.269 -15.125 -15.125 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.12503    8.78569   1.722    0.146
## O2_uM       -0.04187    0.10411  -0.402    0.704
## 
## Residual standard error: 19.85 on 5 degrees of freedom
## Multiple R-squared:  0.03134,    Adjusted R-squared:  -0.1624 
## F-statistic: 0.1618 on 1 and 5 DF,  p-value: 0.7041
#Asv2341
m.norm %>% 
  psmelt() %>% 
  filter(OTU=="Asv2341") %>% 

  lm(Abundance ~ O2_uM, .) %>% 
  summary()
## 
## Call:
## lm(formula = Abundance ~ O2_uM, data = .)
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  9.6487 -2.0548 -1.8821 -2.3513 -2.3513  0.7907 -1.8000 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  2.35128    2.13011   1.104     0.32
## O2_uM       -0.01450    0.02524  -0.575     0.59
## 
## Residual standard error: 4.812 on 5 degrees of freedom
## Multiple R-squared:  0.06193,    Adjusted R-squared:  -0.1257 
## F-statistic: 0.3301 on 1 and 5 DF,  p-value: 0.5905
p.adjust(c(0.5905, 0.5924, 0.0001264, 0.00003552, 0.6734, 0.8889, 0.5905, 0.5905, 0.0000101, 0.8889, 0.8377, 0.7813, 0.6851, 0.09034, 0.5905, 0.8984, 0.799, 0.7065, 0.8889, 0.5905, 0.5905, 0.7762, 0.6649, 0.9432, 0.9432, 0.3544, 0.7041, 0.5905), method="fdr")
##  [1] 0.943200000 0.943200000 0.001179733 0.000497280 0.943200000
##  [6] 0.943200000 0.943200000 0.943200000 0.000282800 0.943200000
## [11] 0.943200000 0.943200000 0.943200000 0.632380000 0.943200000
## [16] 0.943200000 0.943200000 0.943200000 0.943200000 0.943200000
## [21] 0.943200000 0.943200000 0.943200000 0.943200000 0.943200000
## [26] 0.943200000 0.943200000 0.943200000

Plots to go along with these tests.

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>% 
  psmelt() %>% 
  
ggplot() +
  geom_point(aes(x=O2_uM, y=Abundance)) +
  geom_smooth(method='lm', aes(x=O2_uM, y=Abundance)) +
  facet_wrap(~OTU, scales="free_y") +
  labs(title="Example 11: Abundance of ASVs within Family Oceanospirillaceae across oxygen")

m.perc %>% 
  subset_taxa(Family=="D_4__Oceanospirillaceae") %>%
  psmelt() %>% 
  
ggplot() +
  geom_point(aes(x=O2_uM, y=OTU, size=Abundance, color=OTU)) + 
  scale_size_continuous(range = c(0,5)) +
  labs(y="ASV", title="Example 12: Abundance of ASVs within Family Oceanospirillaceae across oxygen")+
    theme(plot.title = element_text(size=15, hjust=.5,vjust=.5,face="plain"),
        axis.text.x = element_text(colour="grey20",size=10,angle=0,hjust=.5,vjust=.5,face="plain"),
        legend.box.background = element_rect(),
        legend.box.margin = margin(6, 6, 6, 6),
        legend.position = "right")

#metadata = data.frame(mothur@sam_data)
m.meta.alpha %>% 
  arrange(Depth_m) %>% 

ggplot(aes(x = O2_uM, y = Depth_m)) +
  geom_point() +
  geom_path(aes(group = 1)) +
  scale_y_reverse(lim=c(200,10)) +
  theme_classic() +
  #options(repr.plot.width=5, repr.plot.height=3) +
  labs(y = "Depth (m)",
       x = "Oxygen (uM)")